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Abstract 


Image processing by partial differential equations (PDEs) has been an ac- 
tive topic in the area of image denoising, which is an important task in 
computer vision, In PDE-based methods for unprocessed image pracess= 
ing, the original image is considered as the initial value for the PDE and th 
solttion af the 


squation is the outcome of the model. Despite the advan- 
tages of using PDEs in image processing, designing and modeling different 
equations for various types af applications have always been a challenging. 
fad interesting problem, In this article, we aim to tackle this problem 
by introducing a fourth-order equation with flexible and trainable coeffi 
Gents, and with the help of an optimal control problem, the coefficients are 
determined; therefore the proposed model adapts itself to each particular 

is performed on the 
noisy test image and the performance of our proposed method is compared 
to other PDE-based models, 


application. At the final stage, the image enhancemes 
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1 Introduction 


Image denoising is an important preliminary step in many computer vision 
and image processing problems that can affect further processing, and has 
been an open research area for a long time [28, 31, 88]. By considering the 
image as a fmction J in R2, the aim of denoising an image is to extract the 
clear image I from the noisy image Jy that has been degraded by the model 
Jy = 1 +n, in which n is a noise fimetion. 

Various approaches for image denoising have been developed, such as 
methods based on nonlocal means filters [4, 18], wavelets [8, 25], Perona and 
Malik (PM) [29], block-matching, and three-dimensional filtering (6, 7, 26] 
sparse representation [46], adaptive image filtering (1, 2), bilateral filtering 
(35, 10], and so on. 

In contrast to the discrete methods used in signal processing in the past, 
which were the basis of digital image processing, in the past decades, con- 
timions methods based on partial differential equations (PDEs) have become 
an effective tool. 

‘The first efforts of using PDEs in computer vision and image processing 
tasks date back to the 1960s [I1, 15]; however, these methods did not attract 
much attention. In the 1980s, the importance of multi-scale descriptions of 
images has been recognized. The importance of scale-space filtering in images 
was introduced by Witkin [41] and then developed by Koenderink [17]. The 
idea of their approach is to generate a family of sealed images I(r,y,f) by 
convolving the original image Jo(x,y) with a Gaussian kernel G(r, y.t) of 
variance t as follows: 


F(x,y,t) = Lul(x,y) 


Glx,y,t) 


‘The time variable # is said to be the seale parameter, and the higher the 
scale gets, it results in an image with reduced resolution; therefore, the noise 
vanishes. 

It was pointed out by Koenderink (17] that, (x,y, ) is the fumdamental 
solution to the heat conduction equation as below, with the original image 
Jo(c,y) a8 the initial condition and the conductance coefficient e, a real num- 
her: 

Ty = c(Tcx + Ty). 


Although the outcome of this equation is a less noisy image, the main 
disadvantage of this approach is the fact that while smoothing, it blurs im- 
portant image features; in other words, the Gaussian scale-space does not 
respect the edges of objects im an image. Therefore, the scaled image loses 
important details after the process 

In 1992, a method introduced by Perona and Malik on anisotropic diffu 
sion [29] drew great attention to PDE methods in image enhancement. In 
fact, they designed the following equation, known as PM second-order PDE, 
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to achieve a good balance between noise removal and retaining the edges of 
objects in the image: 


Ty = e(2,4,t) Lee + Iyy), 

In the above PDB, the conductance coefficient fnnetion, e(1,y,#), is cho- 
sen to decrease smoothing around the edges and reduce noise in other areas. 
Following their model, many equations have been introduced on this basis: 
see for instance (16, 39, 5]. 

Although the PM second-order PDE and its variances were sticcessfully 
applied in image denoising problems; these methods cause blocky effects on 
the processed image. In [43], it is noted that the PM PDB is a second-order 
model, and this feature guarantees its ability to reconstruct images with 
discontinuities but is responsible for the blocky effect [45] 

In recent years, many efforts have been made to resolve this problem in 
the PM method, stich as introducing different. coefficients, different equa- 
tions for processing, raising the order of the equation, and so on. One of the 
innovative efforts to tackle the problems of second-order approaches was in- 
spired by learning-based methods in machine learning {19}. Afterward, some 
second-order learning PDEs (L-PDEs) were introduced that adopt a tech- 
nique called PDE-based optimal control {14] and use training images to learn 
the coefficients in the PDE, and once they are computed, then the PDE is 
obtained and can be applied to test images. 

In learning PDEs proposed in (19, 22, 23], the coefficient functions are 
assumed to minimize a functional subject to some PDE constraints that are 
designed according to the problem. 


In these classes of approaches, the structure of the problem is considered, 
and the learned coefficients will form the PDE according to the training 
images. ‘Thus, the most effort in obtaining a PDE is to prepare some in- 
put/output training image pairs 

One of the other successful strategies in enhancing the performance of 
second-order PDEs for image denoising is to design higher-order equations in 
replace of quadratic ones (13, 21, 24]. In 2000, in a method proposed by You 
and Kaveh, known as YK fourth-order PDB [42], the following equation was 
introduced: 


al 2 

Be TM Ces (19H) 7A), a 
where V2 denotes the Laplacian operator and conductance coefficient (+) 
in YK fourth-order PDE, is a function of the absolute value of the Laplacian 
of the image intensity, that is, 


eo 
oI) = ee @ 
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in which k > 0 is the Laplacian threshold. ‘The YK fourth-order PDE at- 
tempts to remove noise and preserve edges by approximating an observed 
image with a piecewise planar image [42]. 

In recent years, more PDE-based models have been developed, improving, 
the performance of the YK method. In 2016, Zhang and Ye [47] proposed 
the adaptive fourth-order PDE for noise removal 


Ou div h(Vu, Du) 


+ Dal ra) + Haid = h( Wu, DPu)) 


where h( Vu, D®u) is an edge detector chosen by 


in which + > (is a small constant to guarantee h <1, ky and hy are positive 
constants, and jx > 0. and A > 0 are parameters balancing the contribution 
of the three terms in the objective function. 


Later in 2018, Siddig et al. [33] introdnced the following equation: 


=0, (2,1) € Or = (0,7) x 2, 


(x,t) € (0,7) x aD, 
(x,t) € (0,7) x 0, 
ren, 


where a(u) = Foto and Gp(z2) is the Gaussian filter with as its 


parameter, In these models (Siedig’s and Zhang's PDEs), they tend to over- 
come the disadvantages of YK and TV model while keeping the pros of PDE 
algorithms. 

‘The theoretical analysis in (9, 12], shows the advantages of fourth-order 
equations over second-order ones. First, fourth-order diffusion reduces fluctu- 
ation faster than second-order ones. Second, since the second-order PDE will 
evolve toward a piecewise constant approximation in smooth regions, unlike 
second-order PDE, fourth-order PDE will evolve toward and settle down to 
1 piecewise smooth image if the image support is infinite [44] 

In this article, motivated by the previous efforts on solving second-order 
equations’ problems, we propose a fourth-order PDE with trainable coeti- 
cients for image denoising to enhance the performance of the second-order 
L-PDEs. To mention two advantages of this PDE, one should note the ad- 
vantage of using fourth-order PDEs, and also the proposed PDE uses image 
differential invariants up to fourth-order, which can detect more image fea- 
tures compared to those of second-order PDEs. The results obtained from 
our scheme are compared to the state-of-the-art previons models, such as 
second-order learning-based PDEs [22], YK fourth-order PDE [42], and a 
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more recent fourth-order PDEs, Siddig [33], and Zhang's [47], om this topic. 
‘The experimental resnlts show that, in terms of subjective and objective 
measures, the proposed model can ontperform the pre-existing methods, 

‘The rest of this article is organized as follows. Section 2 presents the 
trainable fourth-order PDE with a detailed description of the solution of the 
optimal control problem. Experimental results are presented in Section 3, 
and the article is concluded in Section 4. 


2 Proposed model 


In this section, first, the formulation of the fourth-order PDE in our model 
for noise reduction is presented and then an optimal control problem will be 
introduced to assist in completing the PDE. Finally, the optimality conditions 
will be investigated for the recommended model, which helps us to extract 
enough information to solve the problem. 

‘The framework of this article is based on the fourth-order noise removal 
PDE models, and inspired by second-order learning PDEs mentioned in the 
previons section we will consider flexible trainable coefficients for the pro- 
posed PDB. 

In this work, the typical notations in image processing and optimal control 
literature are used. We denote f as the input image and J as the desired 
output image. Moreover, is an open bounded region of R?, 20 is its 
boundary, the spatial variable (7, y) belongs to 0, ¢ € (0, 7,) is the temporal 
variable, © x (0, Tj) is named Q, and T is its boundary. We then define the 
set 17 as bellow and |p| is the length of string p in which p € 7p 


n= {Ory ry. 22, yy cre, rey. ryy, yyy cern, rery, rryy,cyuy, YUU} 


Changing the viewpoints of an object causes a geometric deformation, 
which can be modeled by a two-dimensional Enclidean, similarity, affine, 
or projective transformation group. Many methods have been designed to 
correctly recognize planar objects, by extracting image invariant features for 
the action of various two-dimensional transformation groups, and differential 
invariants are one of them. We can express image differential invariants as the 
functions of image partial derivatives, and they are widely used to describe 
local image structures. Therefore, to construct our PDE, we used a set of 
differential invariants up to fourth-order. We will use a set of differential 
invariants that are more wildly used. 

‘Thus, the main reason for using these types of partial derivatives to form 
our equation is to create a model that can cover more image features com- 
pared to second-order ones and therefore be more sensitive to edges in the im- 
age. Furthermore, these differential invariants stay unchanged under varions 
two-dimensional transformations like rotation and translation [27]. Consider 
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the fourth-order PDE below: 
{Rima in Q, 

I(zyt)=0, on P, «) 
Ha,y.0) =f. in 2 


in which u(t) = [uo(t) --- ue(t)]” is the coefficient function and diff(-) in- 
cludes the differential parts defined as 


aif (1) = 1, 
diffe) = 1? + 2, 

iffg(1) = Lee Ty — Tey”, 

eiffy(1) = Feely? — Qeylely + Iyyle?, 

diff(1) = Leg Dely — Leyly* — Leyla? + Iyylely 

diffe(I) = Lnsxle® + Bley lyle? + 3Dzyyly Le + Iyyuly’, 

iffy(1) = Leezele! + AL eyyyly Te + 6lexyyle?Ty? + leery lyde® + Iyyyyly!. 


(3) 
‘The initial value for this equation is f, the noisy image, and the ontput, 
I{.) is the desired clear image. To obtain a workable equation, we need to 
determine the coefficient function u(t), which is used to control the evolution 
of J. We ain to attain a flexible system that can train the coefficients; then 
the function u(t) would vary in different problems and be adapted according 
to the given training images. 
As mentioned previously, by using differential invariants in onr model, the 
PDE system (4) is rotationally invariant and the coefficient functions must 
be independent of (r,y) so that (4) is shift-invariant as well, 


Proposition 1. Considering equation (4), the control functions {1}! must 
be independent of (1,9). 


Proof. Considering 
D(I,u) = u(t)"dift(), (6) 


D(I,x,y,t). Therefore, it is enough to prove the 


we then rewrite D(I,w 
independence of D of (x,y). 

Due to the translation invariance of (4), when we change f(x,y) to f(— 
yy), Ua.y) changes to I(x —2!,y—y') and we have 


au 


=y') 


De —2',y—y'),2,y, 1). 


Wi 


‘Then, by the conversion x — 2! =x and y—y! = y, we get 


Al(x.y) 


D(x,y),2+',y + yt). 


OF 


On the other hand, 
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Al(ey) _ 
OF 


DU (x,y), 0,4, t), 


> Dilx+a'y+yt)=DU,2,y,t) for all (2',y/). 


Hence, D is independent of (1, y), which indicates the independence of 
the function u. a 


2.1 Extracting the coefficient function 


Attempting to solve the problem (4), we need to determine the coefficients 
u(t), which is done by an optimal control problem, Therefore we recollect 
some necessary definitions in this field. We begin by assuming I € H¢(0), 
defined to be the closure of the infinitely differentiable fmnetions compactly 
supported in 9 in Hy(Q) 


Definition 1. Considering the optimal control problem 


min J(I,u) 
st T(I,u) =0 7) 
we Ud 
where J : Hj(®) x L4(Q) + Bis called the objective function of the optimal 
control problem and Ua € £°(0) is the admissible set, a nonempty set in 
Rm" 
A vector @ € Una is called an optimal control for (7), if J(I,a) < J(u) 
for each u € Usy and J is called the optimal state associated with i 


Definition 2 ({s6]). Let U and V be normed vector spaces, let Usa © U. 
and let f : Usa + V, and consider some u € Una. Let h EU. fu +th € Usd 
for sufficiently small t > 0, and the limit, 


S(t th) = flu) 
a Sn 


flu, h) = lim: 


exists in V, then 5 f(u,h) is called the directional derivative of F at u in the 
direction h. If the directional derivative 5 (u,h) exists for each A, then the 
map 

6 flu.) sU + Vhs 6f(ush), 
is called the first variation of f at u. In addition, if the first variation con- 
stitutes a bound linear operator, then it is called the Gateanx derivative of 


fi 


Definition 3 ({36]). Let X and ¥ be Banach spaces and let A: X + Y be 
a bounded linear operator. The map 
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AN LY" 4X", A"(f) = fod, 


is called the adjoint of A, For the optimal control problem (7), the adjoint 
equation is of the form 


D,T(I,u)"A=VyJ(y,u), (8) 


and its solution A is the adjoint state. 


Definition 4 ({36]). The fiunetion L(I,u, A) : H(@) x L2(0) x H}(Q) + R 
defined as bellow is called the Lagrangian fimetion for the problem (7) 


[(L,u,d) = du) — ((L.10), A) (9) 
where (-,-) denotes the inner product defined on 3(0). 


Our motivation here is the optimal heating problem, which uses an op 
timal control problem to get the desired output from a heat equation. For 
instance, in a transient heating problem, the goal is to reach a certain tem- 
perature g by controlling the heating elements [30]. Likewise, our objective in 
this article is to approximate the target clean noise-free image T by control- 
ling the coefficient function. Therefore, a PDE-based optimal control problem 
will be presented. As the first step, we need to prepare some input/output 
training images (fc, Ie), where fr is the noisy image and J, is the clear im- 
age. The final output of the equation needs to be close to the ground truth; 
thus the coefficient functions must minimize the following ftmetional: 


I({Le Ha) 


wf, (Iu(Ty) = Ta)"d + 5 ide fe (edt, (10) 


where KC is the number of the training images, [4(7)) is the output image 
determined from (4) at time t = Ty when the initial value is fe, and ay are 
positive weighting parameters '. The first term of the functional J requires 
the final outpnt of our PDE to be elose to the ground truth and the second 
term is for regularization so that this optimal control problem is well-posed 
and counteracts the tendency of the control to become locally bounded as J 
approaches its minimum [36] 
We have the following optimal control problem with PDE constrains: 


Oh — u(t)" diff(U,), in Q, 
min) Hava)s ste} Ie (ayt) =O, on P, (qu) 
A(z,y,0) = fa, in. 


In addition, we consider an inequality constraint 0 < u(t) for the problem 
(11) to make sure that the coefficients are all positive, In the following, the 


Vin this article, we fix as = 10°74 = 0,...,6 
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necessary conditions for this problem (11) will be examined, so that we would 
be able to solve the problem by proper numerical methods. 


Lemma 1, Suppose that J(J,u) = Jy g(u,1)dQ is of elass C*, where g is a 
smooth fimetion, Then considering the function D(I,u) in (6) as 


DJ 
Da = Dill. w) + gL uyLy), 
where Dj and gj are the adjoint operators of D; and gy, respectively, 24 


is the Gateaux derivative of J with respect to the control u, and A is the 
adjoint state of (4) 


Proof. We first define an operator , mapping 1 to the solution of (4), and 
set J = (u). Now consider 


Ulu teu 


where du is the perturbation of u. 
We have 


a2) 
= f (orn fa) + ae. 
4 
(an( 1), 1) + (a, )(8u), Ag 
(97 (u, DQ). da + (97 (u, 1)(1). du)g. 


If we take as the adjoint state of (11), then from the PDE constrains 
of (11), e 


(Gh) Ng = (D(u,¥(4)), dQ (13) 
and 


Now by subtracting (14) with (13 
((W(ut 26u) — Wu) NQ 
= (D(u + .6u, wu + €.5u)) — D(u,v(u)), Ng (15) 
= (€-Dy(utL)(Su) + Dr (u 1) u(u + 2.5u) — ¥(u),d) + o(e). 


In this step, we divide both sides with ©, and let © +0. Therefore 


ma Khoeiniha, Hosseini and Davondi 
(de A) = (Du(u,t)(6u), gq + (Dr(u, 1)(d).A)qQ: (16) 
Integrating by parts implies 
(4, A)a| — (4, Adg = (84, Di (u,1)(M)q + (d,Dj(u, D(A). (17) 
Now. since ) is the adjoint equation of the problem, we have 
(d, Aol3 = 0. (gf (u,1)(1), dq = (—Ar — Dj (u, D(A), d). (as) 
and combining it with (16) yields 
(97(u, DMP), d)q = (—Ac — Dy(u, 1)(A),d) = (Su, Dy (u1)(A))- (19) 
Therefore we arrive at 
DI 


Da 7 9H INT) + Deu DA) 


Now, we investigate the conditions that the optimal vectors a and J must 
satisfy to be then determined, using numerical methods 


‘Theorem 1, Consider the optimal control problem (11), and suppose that 
Una = {0 € La(Q); Ug <u} (20) 


is the admissible control set, where Usa © O C £3(0), O is an open subset 
of L(). In our ease, uy = 0 and J: Hi(2) x La(Q) + R is of class C1, 
If @ € Ugg is an optimal control for (11) in the sense of Definition 1 with 
corresponding state J and adjoint state A, then & = (ii)),cjc, satisfies for 
each 1 € {0,...,6} sin 


tgys where Ly(I,it,) > 0. (21) 


In addition, by introducing the extended Lagrange fimetion 
E(T,u,d, Ha) = L(L,u, A) + (ta ~ ty Ha) + (22) 


and letting 
Ha == max{0, Lu(I,u,d)}, 

the d-tuples (7, 7,\, jiu) © (HG(O), £2(0), H3(O), L2(Q)) satisfies 
K 


ad 
ye l, Agdiff,(Ik)dQ2 = fa, 1 =0,...,6. (23) 
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Proof. First note that, (20) implies that Ugg is convex, We have the following 


optimality system which can be used to determine the optimal control and 
I: 

VaLli,ti,A) = 0, 

VLG, X) = 0, 


VuLlg,i,d) =0, (24) 


(Lu( Lit, A). =a) > 0, for each w € Ua 
Every solution to the general form of optimal control (7) must satisfy this, 
system [36]. 


Slightly rearranging the fourth optimal condition, we get 
(LaF), @) pag) S (Lu(T,,4),0) paca) 


It implies that i is the solution to the optimal problem 


ane (Ea(T, 0,0.) p2¢0) = min 7 Lu(T, i, X)us 5) 


‘The components u, can be varied independently becanse of the form of Usa, 
such that the sum in (25) attains its min if and only if each summand is 
minimal. Therefore, 


LE 4,3) = sn Lay (Ly iy A) (26) 
Now, (21) is a direct result of (26) 
Moreover, according to the definition of L(I.1,., a) in (22), we have 
VaL(L, ty Ay Ha) = VL, i, A), 


ViL(I,u,A, Ha) = ViL(T, i, d), (27) 
Vul(lsu, Asa) = Vu(L,t,\) ~ poe 


‘Therefore, we contimne the proof with the old form of Lagrangian for simplic- 
ity. By Definition 4, the Lagrangian function for the optimal control problem 
(11) is as follows: 


EU) = 


d+ fry, ~ Dile.w))aQ. (28) 
im 


where the multiplier A, is the adjoint state, It can be seen that the first 
optimality condition Vs L(g, i, A) = 0 is the PDE constraint in (11) 


‘To find the adjoint state A, we first perturb D(J,) with respect to I: 
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Dll +.5h.u) — Dll.) 


_ (22 16h) aD 6K) a 
= 6 (tae + at ae) TOO (29) 
=2 Deptt 046 
in which 
— OD) _ yr yy) 2¢@ifEUR)) _ gif) og 
colle) = SR OS a), 


turb Ty in L: 


Loves Tey-+) 


[iG +ebh)auTy) —Ta(e.y)) do 
f 


= 5 [ Coen) ~ Ieemayyan 
i 


+ fala tosh), = DUI + 2:61:)) ~ ((U), — DUe))]4Q) 
é 


= [ten Ty) — In (a,y))8le(x, y, Ty)dQ 


a 
ab 
+ fw (5e(x,y,T))),dQ — [Eo ois) 


dQ +02), 


and 6J, should satisfy 
5h =0 on, 


bh =0 ind 


Due to the boundary and initial conditions of Ig, if we assume that Xx = 0 
on P, then upon integration by parts (with respect to in (J,), and spatial 
variables in 0!?!(6J,)), we have 


Li (lod) = f Can) ~ Halos ihacesn TAM + f Oud yey Tae 


a a 


[9s 5hdQ — [ro sy 2d) ee 


Qe 


bldQ 
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-{ [ve + Te = 1e)6(¢= Ty) = vw), — (ne snag +O(e), 
a 


where d(-) is the Dirac function, Finally, by considering the second optimality 
condition Ly(J,u,A) = 0, our argument yields the following system as the 
adjoint equation for Ae 


2M + T (-1)"\apde)y = 0, in Q, 
are on, (31) 


de(esusTy) = Te—e(Ty), in O, 


for k=1)...,K. 
Now by perturbing u, in L, we get Ly,(,u,\) 


eMiyese) 


v 


pou [fl oresuttoa xf u2(é)de 


Ae((us + .dus)(t) — uj(t) dif f(1)dQ) 


rale 


. rk 
= | (aiu.du,)(t)dt — ( Agdif f (I, )dQ) dude. 
/ [3 [rain 


‘Therefore, we get 


« 
Lu, (Lu,d) = aie, -yf Neti, (L)d, 
ila 


where the function Ay is the solution to (31). Now we get back to the La- 
grangian (22) of our model, we ha 


K 


Ly (Let ea) = 0.0 = S> ft duff U)d = 


a 
rt 


+6, 


e BA = Ly (I,u, A). Therefore we arrive at 


From Lemma 1, we 


Dd > | 
SY Sam Auciff,(T)d0 — pas, 7 =0,- 
Da > [Nelli (Te)aO— 
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In summary, we have the Gateaux derivative of the object function J with 
respect to the control a, 


4 a 
Bag wet Df NH FIL Hass =D 
where the adjoint state ) is the solution to the PDE bellow: 
2M + 3 (-1)'(o,A,), = 0, in Q, 
Ne =O, on, (34) 


Nay Ty) =Te—h(Ty), in 


For more details and a more mathematically precise explanation, the in~ 
terested reader is referred to [36, 20]. Now that we have extracted enough 
information from the optimality conditions, we can use numerical methods 
to determine a for (11). 


Now, we need to find the initial u(f) as the starter of the numerical 


method. At each time step, 24 js expected to be 4-40 so that J, tends 


to the expected output J,. On the other hand, we have 450 = L(J,,u) 
Next, we aim to find u(¢) to minimize 


¥ fut (LU, 0) — td, a-¥ f fet)" u(t) — de(t))"a, 


‘Therefore, the initial u(t) can be obtained by solving the following system: 


P(t)u(t) = d(¢), (35) 


x x 
where we have P(t) = 3° fr pe(t}pn(t)" df and d(t) = 3° frpe(t)de(t) "de, 


in which pg(t) = diff(Iy) and dy(t) = eee 

Finally, with the help of the above u(t) and by considering the Gateaux 
derivative (33), we can the perform Conjugate Gradient (CG) method [32, 34] 
to solve the optimal control problem (11) 


Our PDE framework is summarized in Algorithm 1, Once the coefficients 
are learned, equation (4) is completed and can be applied on test images, 
where the noisy image f is the input as the initial value and the output 
I(T) is the desired denoised image. 


‘Trainable fourth-order pattial diferential equations for image noise removal 249 


‘Algorithm 1 The framework to learn PDEs for image westoration 


Input: Training image pairs (fy,1x), k = 1,...,K , Ty. 
1: Initialize u(t), ¢ € (0, 7)), by solving (35). 
2 while not converged do. 

Compute 22,1 =0,...,6 using (33). 


i 
4: Perform conjugate gradient method using 2 (32, 34). 
5: end while 
6: return Coefficient functions u(t), t € 0, Z)), 
7 Solve (4) using the result from previons step. 


3 Experimental Results 


In this section, the application of our proposed PDE for image denoising is 
demonstrated. ‘The experiments are performed on gray-scale images. We 
compare the results obtained from our PDE with the evolutionary previous 
works, PM second-order PDE [29], YK fourth-order PDE [42], second-order 
L-PDE {22}, Zhang [47], and Siddig’s [33] fourth-order PDE approaches for 
image denoising. For selecting the parameters in PM and YK method, we 
chose the values that produce the most appealing results and the best result is 
chosen as the output of the method. For all experiments, the grid sizes ht = 1 
and Ty = 5 are considered. All our experiments are performed in MATLAB 
R2018b on an Intel(R) Core(TM) 15, 2.40 GHz, 8 GB RAM laptop. 

For a quantitative comparison between the previous methods and our 
scheme, the peak signal-to-noise ratio (PSNR)|44] and structural similarity 
index (SSIM)[37] of the processed images are considered. It should be noted 
that PSNR is most easily defined via the mean squared error (MSE). Given 
‘a noise-free m xn monochrome image J and its noisy approximation f, MSE 
is defined as: 


Loe 


MSE= OT Wa) — 14.3)" 
‘Then PSNR is obtained as follows: 


PSNR 


(36) 
= 20.logyo(max I) — 10.logyo( MSE), 


Here, max I is the maximum possible pixel value of the image, It can 
easily be seen that, the closer an image is to its ground truth, the higher its 
PSNR is. Given similar assumptions, SSIM is defined as 
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SSIM = —y2HE+ CQan 40s) 
Ue ad + Ciley FoF + 


in which [and f are the mean values of I and f, ay and oy are the standard 
deviation of them, 0} and 0? are their variances, and ay, is the variance of J 


and f. For stabilizing the division with the weak denominator, C; = (K1L)* 
and C) = (K2L)* are two variables and L is the dynamic range of pixel 
values, Ay = 0.01 and Ky = 0. Like PSNR, the closer an image is to its 
ground truth, the closer the SSIM will be to one. 

For choosing images in our experiment, first, we used a group of images 
containing 50 clear images and their noisy version, with additive white zero 
mean, Gaussian white noise with the variance of 0.01, in various sizes, cate- 
gorized in ten different classes including animals, sea, personal picture, city 
‘view, nature images, and so on. These images are gathered from SIPI image 
database [40], Berkeley image database [3], and Matlab’s image library. Four 
pictures in each class were used for training the PDE, and one was used as 
the test image and the same test image was used as the initial image for the 
other denoising methods. 

In Figures 1 and 2, the performance of our proposed method and the other 
five denoising models are shown. First, “Lena” picture with size 512 x 512, 
Figure 1(a) is the original clean image, Figure 1(b) is a noisy version of it with 
additive white noise followed by Figures I(c-g) showing the denoised images 
by PM, YK, 2nd L-PDE, Zhang’s’s PDE, and Siddig’s model, respectively, 
and Figure 1(h) is the denoised image with our proposed method. It can 
be seen that compared to the other five methods, our fourth-order PDE 
demonstrates better results in preserving edges while decreasing the noise of 
the image. In Figure 2, the * image with size 850 x 480 is displayed 
followed by Figure 2(b), its noisy version, and on Figures 2(c-g), the enhanced 
form of it by the Pm, YK, 2nd order L-PDE, Zhang's model, and Siddig's 
approach is shown, and in Figure 2(h), the result of our denoising model is 
illustrated. It is vivid that the recommended PDE reduces noise without 
blurring the image features or leaving any speckles. 

The PSNR and SSIM values collected from denoising ten images with 
additive white Gaussian noise with all the mentioned models are presented 
in Table 1. Comparing the values obtained from discussed methods in each 
image shows that our proposed scheme has the best performance among the 
other methods, which is indicated by the higher PSNR and SSIM values. 

In addition, to better investigate the performance of our suggested method 
‘on more complex noises, the second group consisting of 50 images, in five 
separate classes with unknown noise, meaning a combination of zero-mean 
Gaussian white noise, Poisson noise, and the salt & pepper noise (d = 0.1) 
jis added to the images, such as aerial and texture images are chosen. This 
group of images is also gathered from SIPI and Berkeley image databases 
and pictures are taken by a Canon EOS 750D camera. Table 2 represents 
the results of performing all the methods on our five image classes, which 
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can be seen that our snggested model results in better outcomes. In training 
the L-PDEs in [22] and our model for denoising the images with unknown 
noise, we used more training images to learn the coefficients due to the more 
complex noise in each image. In Figures 3 and 4, the noisy test image“Aerial” 
and “Tile” are shown with their denoised versions and it is observed that our 
model is superior in decreasing noise while preserving the image features, 


@ (b) cy 


@) i) c) 


(®) wy 


Figure 1: Comparison of proposed scheme with five methods on “Lena” image with 
auied white Gaussian noise. (a) Original image; (b) Noisy image; (c) PM second-order 
PDE; (4) YK fourth-order PDE: (e) second-order L-PDE; (6) Zhang's fourth-order PDE; 
(4) Sidadig's model; (h) proposed fourtlvorder PDE: 
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(a) (b) © 
(a) (©) 0 
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Figure 2: Comparison of proposed scheme with five methods on *City® image with 
added white Gaussian noise, (a) Original image; (b) Noisy image; (e) PM second-order 
PDE; (d) YK fourth-order PDE; (e) secon-order L-PDE: (0) Zing’ fourth-order PDE: 
(a) Skddig’s model; (h) proposed fourth-order PDE 


In Figure 3, a noisy aerial image with unknown noise is shown and its 


noised forms are Figures 3(b-g). Even with the complex unknown noise, 
the performance of our proposed quadratic model is better than the other 
competitors. In Figure 4, a noisy tile image is presented. As expected from 
the PSNR and SSIM values in Table 2, the recommended learning-based PDE 
still yields better results in comparison to the other contestants. 


Finally, we compare the average PSNR values resulted from all the mod: 
els, in Figure 5. Comparing the average values obtained from each method 
on both image groups with Gaussian and unknown noise, we can see that 
the average values gained from our proposed PDE are larger than the other 


approaches on both types of noises, which indicates that our scheme ontper 
forms other models in both image categories. 
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‘Table 1: PSNR and SSIM results for ten test images with 
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additive white Gaussian 


YK 2nd PDE Zhang Siddig Proposed 
Tena TSG FSIS IT TSsl INTIS 80.4875 
0.8031 0.8093 0.8069 0.8421 0.9184 
City 20.3853 20.7052 24.2705 «28.2597 29.6089 29.1105 
0.6008 0.7288 0.8041 0.8172 0.8133 0.8563 
Text 18.0646 21.9834 23.9050 28.0671 28.5636 31.2598 
0.7093 0.7313 0.8527 0.8530 0.8600 0.8697 
Woman 19.8594 20.7346 0504 29.4682 30.1833 30.4666 
0.6551 0.6995 0.8211 0.873 0.8237 0.9030 
Eiffel 19.1024 21.4389 24.1924 28.0046. 30.1453 80.2713 
0.6694 0.6890 0.7806 0.8074 0.8497 0.8901 
Sea 2.2217 21.0962 23.9878 27.2673 7.H078 1.0694 
0.7589 0.7854 0.8030 0.8864 0.7800, 0.8751 
Nighteity 19.8372 21.9698 25.2569 28.5490 28.0172 29.3441 
0.7048 0.7589 0.8497 0.7573 0.8628 0.8638 
Nature 20.2255 22 24.6548 29.7069 29.6567 30.4155 
0.6508 0.7045 0.7886 0.8349 0.9153. 
Flowers 19.3536 20.9230 25.2612 28.9910 29.0531 29.9513 
0.6461 0.6934 O.7914—0.8199 0.8276 0.8785 
Trees 20.3819 22.1760 25.8049 28.6375 28.8983 30.5744 
0.6835 0.7412 __0.7364___0.8536 0.8619 0.9143 
‘Table 2: PSNR and SSIM results for five test images with real unknown noise 
PM YK and I-PDE Zhang Siddig Proposed 
Trick To116s 172001 18.518 10.9679 21.0701 22.5285 
0.7001 0.7058 0.7413 0.7588 0.8103 0.8526 
Aerial 15.9888 18.1222 18.9908 19.5223. 20.9317 22.9611. 
0.6846 0.7140 0.7851 0.7987 0.8159 0.8794 
Personal 16.6055 18.3664 19.0258 20.6046 21.6548 23.7669 
0.6122 0.6988 0.7423 0.7642 0.7941 0.8682 
Texture 15.4097 17.0267 18.9516 19.5583 20.9513 22.0304 
0.7085 0.7631 -——«0.7970—O.8155 0.8289 0.8437. 
Tile 16.2849 17.9806 18.9823 20.7140 21.1675 23.6656 
0.6827 0.7125 O.7681__—O.S155 8468 0.8549. 
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Figure 
real unknown neise. (a) Noisy image: (b) PM second-order PD! 
PDE; (d) seconorder L-PDE; (e) Zhang's fourtieonder PDE: (f) 
propose fourth-order PDE 


{; Comparison of proposed scheme with five methods on “Aerial” image with 
(©) YK fourth-order 
iddig’s model; (x) 


4 Conclusion 


In this article, we proposed a fourth-order PDE for image denoising with 
trainable coefficients asserted by an optimal control problem. The experi- 
mental results show improvements in performance in comparison to previons 
approaches in the field of PDE-based methods for image denoising. Compared 
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to previous methods, our scheme results in better outcomes, clearer image 
with higher PSNR and SSIM. For future studies, we want to strengthen our 
approach in the following areas. First, design a more general equation to 
adapt better to different problems, more precisely, having an optimal con- 
trol problem with controls depending on spatial variables as well. Second, 
to introduce a region-based PDE with trainable coefficients, which leads to 
having two controls for each region. Finally, finding a more suitable mmeri- 
cal method for solving the optimal control problem to have a faster and more 
accurate numerical result. 


(a) (b) © 


@ 


Figure 4: Comparison of proposed scheme with Sve methods om “Tile” image with 
real unknown noise. (a) Noisy image: (b) PM second-order PDE; (c) YE fourth-order 
PDE; (d) second-order L-PDE; (e) Zhang's fourth-order PDE; (f) Siddis's model; (x) 
proposed fourtlvorder PDE, 
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Figure 5: Comparison of average PS 
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